library(Seurat)
library(tictoc)
library(dplyr)
library(ggplot2)
library(pheatmap)
library(cowplot)
load("Colon_0.6_0722.RDS")
levels(immune.combined)
In this step, we cluster the cells within each cell type again to separate the doublets as outliers. In Some clusters, like cluster 0 and 1, doublets does not form a unique isolated cluster while others do. We remove those doublet minor populations to enhance the doublet removal performance.
tic()
celltype='2.CD4+CTLA4+ T'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()
options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
DimPlot(sub, group.by="DF")+ggtitle(celltype),
FeaturePlot(sub, features="doublet"),
DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
ncol=3
)
doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='7'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"
tic()
celltype='3.CD4+ Teff'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()
options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
DimPlot(sub, group.by="DF")+ggtitle(celltype),
FeaturePlot(sub, features="doublet"),
DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
ncol=3
)
Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 1165 Number of edges: 55853 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.5346 Number of communities: 8 Elapsed time: 0 seconds 4.43 sec elapsed
doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='7'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"
tic()
celltype='4.Myeloid phagocytes'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()
options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
DimPlot(sub, group.by="DF")+ggtitle(celltype),
FeaturePlot(sub, features="doublet"),
DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
ncol=3
)
Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 1029 Number of edges: 39231 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.6309 Number of communities: 6 Elapsed time: 0 seconds 4.022 sec elapsed
doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='4'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"
tic()
celltype='5.Ly2Z+ Myeloid'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()
options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
DimPlot(sub, group.by="DF")+ggtitle(celltype),
FeaturePlot(sub, features="doublet"),
DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
ncol=3
)
Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 756 Number of edges: 27536 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.6967 Number of communities: 6 Elapsed time: 0 seconds 3.17 sec elapsed
doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='5'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"
tic()
celltype='6.IgA+ Plasma B'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()
options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
DimPlot(sub, group.by="DF")+ggtitle(celltype),
FeaturePlot(sub, features="doublet"),
DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
ncol=3
)
Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 755 Number of edges: 31786 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.6785 Number of communities: 7 Elapsed time: 0 seconds 3.235 sec elapsed
doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='2'|RNA_snn_res.0.8=='6'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"
tic()
celltype='14.B plasmablast'
sub=subset(immune.combined, celltype0627==celltype)
DefaultAssay(sub)="RNA"
sub=FindVariableFeatures(sub, selection.method = "vst", nfeatures = 1500)
sub=ScaleData(sub, features = rownames(sub), verbose = F)
sub=RunPCA(sub, features = VariableFeatures(object = sub),verbose = F)
sub=RunUMAP(sub, min.dist = 0.2, n.neighbors = 10, dims = 1:30, verbose = F)
sub=FindNeighbors(sub, dims = 1:20)
sub=FindClusters(sub, resolution = 0.8)
toc()
options(repr.plot.width=20, repr.plot.height=6)
plot_grid(
DimPlot(sub, group.by="DF")+ggtitle(celltype),
FeaturePlot(sub, features="doublet"),
DimPlot(sub, group.by="RNA_snn_res.0.8", label=T, repel=T, label.size=10),
ncol=3
)
Computing nearest neighbor graph Computing SNN
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck Number of nodes: 240 Number of edges: 9463 Running Louvain algorithm... Maximum modularity in 10 random starts: 0.5096 Number of communities: 4 Elapsed time: 0 seconds 1.916 sec elapsed
doublet.barcodes=Cells(subset(sub, RNA_snn_res.0.8=='3'))
immune.combined@meta.data [doublet.barcodes, "DF"]="Doublet"
options(repr.plot.width=10, repr.plot.height=8)
DimPlot(immune.combined, group.by="DF")
# Rename cluster id
Idents(immune.combined)<-"SCT_snn_res.0.6"
new.cluster.ids <- c(
"0.IgD+ mature B",
"1.IgM+IgD+ mature B",
"2.CD4+CTLA4+ T",
"3.CD4+ Teff",
"4.Myeloid phagocytes",
"5.Ly2Z+ Myeloid",
"6.IgA+ Plasma B",
"7.Gata3+ ILC2,ILC3",
"8.RORga+ ILC2,ILC3",
"9.Foxp3+ Treg",
"10.CD8+ T, NKT",
"11.IgD+ mature B",
"12.Neu",
"13.Proliferating B",
"14.B plasmablast",
"15.CX3CR1+ Mac",
"16.Epithelial cells",
"17.Endothelial cells",
"18.Fibroblasts")
names(new.cluster.ids) <- levels(immune.combined)
immune.combined <- RenameIdents(immune.combined, new.cluster.ids)
immune.combined[["celltype0627"]]<-Idents(immune.combined)
options(repr.plot.width=10, repr.plot.height=8)
DimPlot(immune.combined, group.by="celltype0627")
# doublet-free
seu=subset(immune.combined, DF=="Singlet")
options(repr.plot.width=10, repr.plot.height=8)
DimPlot(seu, group.by="celltype0627")
colnames(immune.combined[[]])
save(immune.combined, file="Colon_0.6_0730.RDS", compress=T)